This is the analysis of two runs of Acanthophora spicifera salinity and nutrient experiments conducted on the lanai in St. John 616 in November and December 2023. These experiments incorporated two salinity/nutrient levels and two temperature levels.
Packages loaded:
library(lme4)
library(lmerTest)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(strengejacke)
## # Attaching packages (red = needs update)
## ⚠ ggeffects 1.3.4 ✔ sjlabelled 1.2.0
## ⚠ sjstats 0.18.2 ✔ esc 0.5.1
##
## Warnings or errors in CRAN checks for package(s) 'ggeffects'.
## Update packages in red with 'sj_update()'.
#library(mmtable2) no longer available
library(gt)
library(purrr)
library(stringr)
library(tidyr)
library(piecewiseSEM)
library(easystats)
library(magrittr)
Open weight dataset and make columns for growth rate from initial and final weights
acan_growth <- read.csv("../input/acanthophora_growth_r1.csv")
Make a new column for weight change (difference final from initial)
acan_growth$growth_rate_percent <- (acan_growth$final_weight - acan_growth$initial_weight) / acan_growth$initial_weight * 100
Make a new column for daily growth rate from 8 day study (steady growth rate assumed rather than exponential)
acan_growth$steady_growth_daily <- acan_growth$growth_rate_percent / 8
Assigns temperature as a factor
acan_growth$temperature <- as.factor(acan_growth$temperature)
Assigns treatment as characters from integers then to factors
acan_growth$treatment <- as.factor(as.character(acan_growth$treatment))
Assign run as a factor
acan_growth$run <- as.factor(acan_growth$run)
Assign plant ID as a factor
acan_growth$plant_ID <- as.factor(acan_growth$plant_ID)
Assign RLC order as a factor
acan_growth$rlc_order <- as.factor(acan_growth$rlc_order)
Create subsets for the plots and remove one datapoint that is strongly skewing data (overall result does not change)
acan_growth <- subset(acan_growth, secondary_apices != "200")
acan_growth$treatment_graph[acan_growth$treatment == 0] <- "1) 35ppt/0.5umol"
acan_growth$treatment_graph[acan_growth$treatment == 1] <- "2) 28ppt/14umol"
acan_growth$temperature_graph[acan_growth$temperature == 24] <- "24°C"
acan_growth$temperature_graph[acan_growth$temperature == 28] <- "28°C"
Make a histogram of the growth rate data
hist(acan_growth$growth_rate_percent, main = paste("Acanthophora spicifera 9-Day Growth (%)"), col = "maroon", labels = TRUE)
#or
acan_growth %>% ggplot(aes(growth_rate_percent)) +
geom_histogram(binwidth=5, fill = "#990066", color = "black", linewidth = 0.25, alpha = 0.85) +
theme_bw()
##Run model without RLC_order as this has little effect, add initial number of apices
growth_model_acan <- lmer(formula = growth_rate_percent ~ treatment + temperature +
(1 | plant_ID) + (1 | run) + (1 | initial_axes), data = acan_growth, REML = TRUE)
Plot histograms and performance checks
hist(resid(growth_model_acan))
plot(resid(growth_model_acan) ~ fitted(growth_model_acan))
qqnorm(resid(growth_model_acan))
qqline(resid(growth_model_acan))
performance::check_model(growth_model_acan)
Get R2 and summarize
rsquared(growth_model_acan)
## Response family link method Marginal Conditional
## 1 growth_rate_percent gaussian identity none 0.008527329 0.7167619
summary(growth_model_acan)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: growth_rate_percent ~ treatment + temperature + (1 | plant_ID) +
## (1 | run) + (1 | initial_axes)
## Data: acan_growth
##
## REML criterion at convergence: 772.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.15114 -0.44427 0.05451 0.46396 1.98607
##
## Random effects:
## Groups Name Variance Std.Dev.
## plant_ID (Intercept) 148.58 12.190
## initial_axes (Intercept) 49.05 7.003
## run (Intercept) 67.40 8.210
## Residual 105.99 10.295
## Number of obs: 95, groups: plant_ID, 48; initial_axes, 8; run, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.864 7.179 1.765 -0.399 0.733
## treatment1 1.229 2.141 44.627 0.574 0.569
## temperature28 -3.348 4.155 38.429 -0.806 0.425
##
## Correlation of Fixed Effects:
## (Intr) trtmn1
## treatment1 -0.141
## temperatr28 -0.288 -0.007
View random effects levels and make a table of the model output
ranef(growth_model_acan)
## $plant_ID
## (Intercept)
## 1 9.53417840
## 2 11.50686729
## 3 3.87228371
## 4 7.48459427
## 5 -8.96283234
## 6 20.49190957
## 7 -15.56038211
## 8 -12.17546306
## 9 10.63486896
## 10 15.98775355
## 11 6.44751682
## 12 -0.91797790
## 13 4.67935909
## 14 6.86528363
## 15 -11.18300107
## 16 6.05450444
## 17 -10.95411555
## 18 13.84194148
## 19 -18.33904191
## 20 -7.22510672
## 21 8.01159244
## 22 -18.95843963
## 23 -7.19580500
## 24 -1.89213433
## 25 -9.34308528
## 26 0.86169829
## 27 11.39784645
## 28 -11.62101683
## 29 0.73408357
## 30 0.08916723
## 31 -10.19593676
## 32 10.63710508
## 33 6.49379722
## 34 12.18624772
## 35 -0.15971334
## 36 -9.20885772
## 37 10.61902370
## 38 3.56617631
## 39 7.37701722
## 40 -4.37985585
## 41 -9.74270835
## 42 -6.44915512
## 43 -5.04609703
## 44 -21.24271841
## 45 -1.29978106
## 46 -1.83057799
## 47 2.66228410
## 48 11.84670284
##
## $initial_axes
## (Intercept)
## 3 3.0315165
## 4 0.7459840
## 5 2.3403631
## 6 0.7300085
## 7 7.3280180
## 8 -11.2173697
## 10 0.2887361
## 11 -3.2472564
##
## $run
## (Intercept)
## 1 5.465689
## 2 -5.465689
##
## with conditional variances for "plant_ID" "initial_axes" "run"
tab_model(growth_model_acan, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
| growth_rate_percent | ||||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | df |
| (Intercept) | -2.86 | 7.18 | -17.13 – 11.40 | -0.40 | 0.691 | 88.00 |
| treatment [1] | 1.23 | 2.14 | -3.03 – 5.48 | 0.57 | 0.567 | 88.00 |
| temperature [28] | -3.35 | 4.15 | -11.60 – 4.91 | -0.81 | 0.423 | 88.00 |
| Random Effects | ||||||
| σ2 | 105.99 | |||||
| τ00 plant_ID | 148.58 | |||||
| τ00 initial_axes | 49.05 | |||||
| τ00 run | 67.40 | |||||
| ICC | 0.71 | |||||
| N plant_ID | 48 | |||||
| N run | 2 | |||||
| N initial_axes | 8 | |||||
| Observations | 95 | |||||
| Marginal R2 / Conditional R2 | 0.009 / 0.717 | |||||
Plot the model effects
plot(allEffects(growth_model_acan))
Construct null model to perform likelihood ratio test REML must be FALSE
acan_growth_treatment_null <- lmer(formula = growth_rate_percent ~ temperature + (1 | plant_ID) + (1 | run), data = acan_growth, REML = FALSE)
acan_growth_model2 <- lmer(formula = growth_rate_percent ~ treatment + temperature + (1 | plant_ID) + (1 | run), data = acan_growth, REML = FALSE)
anova(acan_growth_treatment_null, acan_growth_model2)
## Data: acan_growth
## Models:
## acan_growth_treatment_null: growth_rate_percent ~ temperature + (1 | plant_ID) + (1 | run)
## acan_growth_model2: growth_rate_percent ~ treatment + temperature + (1 | plant_ID) + (1 | run)
## npar AIC BIC logLik deviance Chisq Df
## acan_growth_treatment_null 5 798.17 810.94 -394.08 788.17
## acan_growth_model2 6 799.41 814.74 -393.71 787.41 0.7559 1
## Pr(>Chisq)
## acan_growth_treatment_null
## acan_growth_model2 0.3846
acan_growth_temperature_null <- lmer(formula = growth_rate_percent ~ treatment + (1 | plant_ID) + (1 | run), data = acan_growth, REML = FALSE)
acan_growth_model3 <- lmer(formula = growth_rate_percent ~ treatment + temperature + (1 | plant_ID) + (1 | run), data = acan_growth, REML = FALSE)
anova(acan_growth_temperature_null, acan_growth_model3)
## Data: acan_growth
## Models:
## acan_growth_temperature_null: growth_rate_percent ~ treatment + (1 | plant_ID) + (1 | run)
## acan_growth_model3: growth_rate_percent ~ treatment + temperature + (1 | plant_ID) + (1 | run)
## npar AIC BIC logLik deviance Chisq Df
## acan_growth_temperature_null 5 798.02 810.79 -394.01 788.02
## acan_growth_model3 6 799.41 814.74 -393.71 787.41 0.6119 1
## Pr(>Chisq)
## acan_growth_temperature_null
## acan_growth_model3 0.4341
Simple boxplots for treatment and temperature
acan_growth %>% ggplot(aes(treatment_graph, growth_rate_percent)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.75, size = 5, aes(color = temperature), position = "jitter", show.legend = TRUE) +
labs(x="treatment", y= "9-Day Growth (%)", title= "A", subtitle = "Acanthophora spicifera") +
scale_x_discrete(labels = c("T0) 35 ppt/0.5 μmol N", "T1) 28 ppt/14 μmol N")) +
ylim(-60, 60) + stat_mean() +
scale_color_manual(values = c("azure4", "darkgoldenrod1")) +
geom_hline(yintercept=0, color = "red", linewidth = 0.5, alpha = 0.5) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
acan_growth %>% ggplot(aes(temperature_graph, growth_rate_percent)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.75, size = 5, aes(color = treatment), position = "jitter", show.legend = TRUE) +
labs(x="Temperature", y= "9-day Growth (%)", title= "B", subtitle = "Acanthophora spicifera") +
scale_x_discrete(labels = c("24°C", "28°C")) +
ylim(-60, 50) + stat_mean() +
geom_hline(yintercept=0, color = "red", linewidth = 0.5, alpha = 0.5) +
scale_color_manual(values = c("azure3", "darkgoldenrod")) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
Summarize the means
acan_growth %>% group_by(treatment) %>% summarise_at(vars(growth_rate_percent), list(mean = mean))
## # A tibble: 2 × 2
## treatment mean
## <fct> <dbl>
## 1 0 -2.96
## 2 1 -1.33
acan_growth %>% group_by(temperature) %>% summarise_at(vars(growth_rate_percent), list(mean = mean))
## # A tibble: 2 × 2
## temperature mean
## <fct> <dbl>
## 1 24 -0.543
## 2 28 -3.73
##RUN MODEL FOR SECONDARY APICES
First plot histogram
acan_growth %>% ggplot(aes(secondary_apices)) +
geom_histogram(binwidth=5, fill = "#990066", color = "black", linewidth = 0.25, alpha = 0.85) +
theme_bw()
Run model for secondary apices
apices_model_acan <- lmer(formula = secondary_apices ~ treatment + temperature +
(1 | plant_ID) + (1 | run) + (1 | rlc_order), data = acan_growth, REML = TRUE)
Check performance of model
hist(resid(apices_model_acan))
plot(resid(apices_model_acan) ~ fitted(apices_model_acan))
qqnorm(resid(apices_model_acan))
qqline(resid(apices_model_acan))
performance::check_model(apices_model_acan)
R squared
rsquared(apices_model_acan)
## Response family link method Marginal Conditional
## 1 secondary_apices gaussian identity none 0.2748104 0.6747653
summary(apices_model_acan)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: secondary_apices ~ treatment + temperature + (1 | plant_ID) +
## (1 | run) + (1 | rlc_order)
## Data: acan_growth
##
## REML criterion at convergence: 777
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.94937 -0.57432 0.00079 0.45150 2.21087
##
## Random effects:
## Groups Name Variance Std.Dev.
## plant_ID (Intercept) 52.30 7.232
## rlc_order (Intercept) 76.24 8.732
## run (Intercept) 49.05 7.003
## Residual 144.41 12.017
## Number of obs: 95, groups: plant_ID, 48; rlc_order, 24; run, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 40.566 6.385 2.219 6.353 0.018356 *
## treatment1 21.965 4.338 12.841 5.063 0.000226 ***
## temperature28 -1.007 4.814 22.713 -0.209 0.836247
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn1
## treatment1 -0.336
## temperatr28 -0.377 -0.005
View random effects levels, print table and plot all effects
ranef(apices_model_acan)
## $plant_ID
## (Intercept)
## 1 2.1006896
## 2 -0.4197587
## 3 0.1660575
## 4 0.1660575
## 5 0.4998379
## 6 -3.7009093
## 7 1.5579201
## 8 -3.2729392
## 9 -2.6842007
## 10 3.4068828
## 11 4.1584404
## 12 4.7885525
## 13 4.1142738
## 14 -1.6234315
## 15 -3.4195436
## 16 -0.8990953
## 17 -3.7248190
## 18 2.9963766
## 19 1.2410021
## 20 2.0811516
## 21 -1.3884242
## 22 -6.8493956
## 23 -2.1572453
## 24 -2.1572453
## 25 -5.8447058
## 26 -5.2145937
## 27 -5.4689270
## 28 -5.4689270
## 29 0.7458996
## 30 5.1566842
## 31 -5.9674006
## 32 0.7537950
## 33 2.3927203
## 34 1.7626082
## 35 3.7743900
## 36 -10.2981132
## 37 11.2248719
## 38 0.7230038
## 39 2.2874895
## 40 3.1276389
## 41 5.7628866
## 42 3.2424382
## 43 -6.0742812
## 44 7.1580725
## 45 6.6290198
## 46 -3.2427362
## 47 -3.1714078
## 48 1.0293394
##
## $rlc_order
## (Intercept)
## 1 -2.5208706
## 2 -5.7058278
## 3 -11.1497242
## 4 -9.7538717
## 5 -0.8423165
## 6 -2.9409828
## 7 4.7802381
## 8 -7.1586863
## 9 7.9545288
## 10 -0.9602800
## 11 -0.8439845
## 12 4.4926149
## 13 6.6672924
## 14 5.8770567
## 15 14.3796372
## 16 -4.2787335
## 17 2.6274081
## 18 -2.3146177
## 19 9.4375799
## 20 8.7370447
## 21 -2.1057634
## 22 -5.7356272
## 23 -4.9661900
## 24 -3.6759247
##
## $run
## (Intercept)
## 1 -4.707122
## 2 4.707122
##
## with conditional variances for "plant_ID" "rlc_order" "run"
tab_model(apices_model_acan, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
| secondary_apices | ||||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | df |
| (Intercept) | 40.57 | 6.39 | 27.88 – 53.25 | 6.35 | <0.001 | 88.00 |
| treatment [1] | 21.96 | 4.34 | 13.34 – 30.59 | 5.06 | <0.001 | 88.00 |
| temperature [28] | -1.01 | 4.81 | -10.57 – 8.56 | -0.21 | 0.835 | 88.00 |
| Random Effects | ||||||
| σ2 | 144.41 | |||||
| τ00 plant_ID | 52.30 | |||||
| τ00 rlc_order | 76.24 | |||||
| τ00 run | 49.05 | |||||
| ICC | 0.55 | |||||
| N plant_ID | 48 | |||||
| N run | 2 | |||||
| N rlc_order | 24 | |||||
| Observations | 95 | |||||
| Marginal R2 / Conditional R2 | 0.275 / 0.675 | |||||
plot(allEffects(apices_model_acan))
Construct null model to perform likelihood ratio test REML must be FALSE
acan_apices_treatment_null <- lmer(formula = secondary_apices~ temperature + (1 | plant_ID) + (1 | run) + (1 | rlc_order), data = acan_growth, REML = FALSE)
acan_apices_model2 <- lmer(formula = secondary_apices ~ treatment + temperature + (1 | plant_ID) + (1 | run) + (1 | rlc_order), data = acan_growth, REML = FALSE)
anova(acan_apices_treatment_null, acan_apices_model2)
## Data: acan_growth
## Models:
## acan_apices_treatment_null: secondary_apices ~ temperature + (1 | plant_ID) + (1 | run) + (1 | rlc_order)
## acan_apices_model2: secondary_apices ~ treatment + temperature + (1 | plant_ID) + (1 | run) + (1 | rlc_order)
## npar AIC BIC logLik deviance Chisq Df
## acan_apices_treatment_null 6 821.05 836.37 -404.52 809.05
## acan_apices_model2 7 805.52 823.40 -395.76 791.52 17.525 1
## Pr(>Chisq)
## acan_apices_treatment_null
## acan_apices_model2 2.835e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
acan_apices_temperature_null <- lmer(formula = secondary_apices ~ treatment + (1 | plant_ID) + (1 | run) + (1 | rlc_order), data = acan_growth, REML = FALSE)
acan_apices_model3 <- lmer(formula = secondary_apices ~ treatment + temperature + (1 | plant_ID) + (1 | run) + (1 | rlc_order), data = acan_growth, REML = FALSE)
anova(acan_apices_temperature_null, acan_apices_model3)
## Data: acan_growth
## Models:
## acan_apices_temperature_null: secondary_apices ~ treatment + (1 | plant_ID) + (1 | run) + (1 | rlc_order)
## acan_apices_model3: secondary_apices ~ treatment + temperature + (1 | plant_ID) + (1 | run) + (1 | rlc_order)
## npar AIC BIC logLik deviance Chisq Df
## acan_apices_temperature_null 6 803.57 818.89 -395.78 791.57
## acan_apices_model3 7 805.52 823.40 -395.76 791.52 0.0469 1
## Pr(>Chisq)
## acan_apices_temperature_null
## acan_apices_model3 0.8286
Plots
acan_growth %>% ggplot(aes(treatment_graph, secondary_apices)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.75, size = 5, aes(color = temperature), position = "jitter", show.legend = TRUE) +
labs(x="treatment", y= "Number of Secondary Apices", title= "A", subtitle = "Acanthophora spicifera") +
scale_x_discrete(labels = c("T0) 35 ppt/0.5 μmol N", "T1) 28 ppt/14 μmol N")) +
ylim(-1, 130) + stat_mean() +
scale_color_manual(values = c("azure4", "darkgoldenrod1")) +
geom_hline(yintercept=0, color = "red", linewidth = 0.5, alpha = 0.5) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
acan_growth %>% ggplot(aes(temperature_graph, secondary_apices)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.75, size = 5, aes(color = treatment), position = "jitter", show.legend = TRUE) +
labs(x="Temperature", y= "Number of Secondary Apices", title= "B", subtitle = "Acanthophora spicifera") +
scale_x_discrete(labels = c("24°C", "28°C")) +
ylim(-1, 130) + stat_mean() +
geom_hline(yintercept=0, color = "red", linewidth = 0.5, alpha = 0.5) +
scale_color_manual(values = c("azure3", "darkgoldenrod")) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
Summarize the means
acan_growth %>% group_by(treatment) %>% summarise_at(vars(secondary_apices), list(mean = mean))
## # A tibble: 2 × 2
## treatment mean
## <fct> <dbl>
## 1 0 40.1
## 2 1 61.7
acan_growth %>% group_by(temperature) %>% summarise_at(vars(secondary_apices), list(mean = mean))
## # A tibble: 2 × 2
## temperature mean
## <fct> <dbl>
## 1 24 51.0
## 2 28 50.5